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Abstract 

We study the evolution of domain wall networks appearing after phase transitions in the early 
Universe. They exhibit interesting dynamical scaling behaviour which is not yet well understood, 
and are also simple models for the more phenomenologically acceptable string networks. We 
have run numerical simulations in two- and three-dimensional lattices of sizes up to 4096 3 . The 
theoretically predicted scaling solution for the wall area density A oc 1/t is supported by the 
simulation results, while no evidence of a logarithmic correction reported in previous studies could 
be found. The energy loss mechanism appears to be direct radiation, rather than the formation 
and collapse of closed loops or spheres. We discuss the implications for the evolution of string 
networks. 
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I. INTRODUCTION 



The idea that the symmetries in nature are not respected by the vacuum plays a crucial 
role in the unification of forces. Moreover, any broken symmetries that can be identified 
today were likely to have been restored at high temperatures ]I], 0, |3|. These two facts 
suggest that the Universe underwent phase transitions in its early history. It was realised 
that phase transitions may leave the Universe filled with topological defects 0, ^ , which if 
massive enough could be observed through their density fluctuations || [7| (see also || || 
for reviews). 

Probably the most interesting and important property of defect networks is that they 
seem to exhibit dynamic scaling. This means that they quickly lose memory of their initial 
conditions and evolve towards configurations which can be characterised by a single length 
scale (or perhaps a few JlO|) £. This length scale is thought to increases with time with a 
universal exponent. In a relativistic field theory one can argue from dimensional analysis 
that C,(t) ~ tf(Mt), where M is the mass scale of the defect. The large-scale dynamics 
of defects are independent of M, and so the dynamical scaling exponent can be naively 
estimated as 1. 

This behaviour has been checked by numerically simulating classical field theories for 



domain walls |TT], [T2|, O], gauge strings |14], [15|, [IBL global strings |T7], ITS], global monopoles 



I9|, EUl 1211 and textures I9|, E22l E3I. All the simulations are consistent with the linear 



scaling law over the range of the simulations, but do allow other behaviours: in particular, 
Press, Ryden and Spergel suggested that the results for domain walls would be better fitted 
by£~t/ln(t). 

Since the original naive scaling arguments were put forward, a more quantitative approach 
to the dynamics of domain wall networks has been proposed by one of the authors |25j, |26], |27| , 
which predicts not only the linear scaling law for domain walls, but also the amplitude of 
the relation. Clearly, a logarithmic correction would be a problem for the approach. The 
main purpose of this paper is to give a more accurate numerical determination of the scaling 
law (in 2 and 3 dimensions) and to determine the amplitude, to check the accuracy of the 
theoretical predictions in [^, p7fl . 

Our results for the scaling exponents and amplitude can be found in Tables |, [TT| They 
are consistent with the linear scaling law for £, but the numerically determined amplitudes 
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are higher than the theoretical predictions. A detailed comparison can be found in 



II. DYNAMICS OF DOMAIN WALLS 

Domain walls occur in field theories whose manifold of minimum energy states is topo- 
logically disconnected (see e.g. ||). The canonical example in relativistic field theory is a 
theory of a single scalar field 4>(x), with action 



S 



rf 4 xv^Q^0-V(0)), (1) 



where the potential V is a function with more than one minimum, which we take to be the 
renormalisable form 

V{<P) = \{^-^f. (2) 

In the cosmological context, the metric is taken to have the Friedmann-Robertson- Walker 
form 

where rj^ = diag(l, —1, —1, —1) is the Minkowksi space metric and t is conformal time. The 
field can be conformally rescaled R(t)<f) — > 0, giving an Euler-Lagrange equation 

<j> + 2^0 - V 2 + A0 (f - fx 2 R 2 ) = 0. (3) 

In the broken phase, (/i 2 > 0) there are domain wall solutions in which the field changes 
vacuum over a distance of order M _1 , passing through zero. The solution for an infinite 
static planar wall is well known || and with an appropriate choice of coordinates can be 
written 

= //tanh(Mz). (4) 

At high temperature the mass parameter receives thermal corrections from the fluctuations 
in fields to which is coupled ||1], [|: 

AT) = £ ~ cT\ 

where c is a model-dependent numerical factor. For the pure scalar field theory c = 1/6. 
Thus as the Universe cools the field undergoes a phase transition from (0) = to (0) = 
±/x(T). As this transition happens at a finite rate, the field cannot select the same minimum 



everywhere at the same time, and the Universe divides into domains in which takes either 
positive or negative values. By continuity of the field, these domains must be separated by 
domain walls || in which the field approximates the configuration given by Eq. f| in the 
transverse direction. The initial size of the domains £ is controlled by the rate at which the 
transition occurs and how strongly damped the field is 



The domain walls are mostly in the form of one infinite boundary separating percolating 
clusters of the two vacua |L2], [L3|, [17|] . The subsequent evolution of the field is controlled by 
the dynamics of this infinite domain wall. The wall has a tension of order M 3 , and tries to 
straighten out and lose energy, and in finite volume eventually one or other vacuum will take 
over the whole space and the field thereby reaches equilibrium. One way of quantifying the 
approach to equilibrium is to measure the area density of the domain wall A, or equivalently 
the curvature scale of the wall £ = 1/A, where j25] 



A=W)|V0|). (5) 

At late times, this quantity can only depend on time t and the mass scale M. The fact 
that the wall obeys a Nambu-Goto equation independently of M |J indicates that A — a/t 
purely on dimensional grounds, where a is a constant amplitude. This can be called the naive 
or "classical" scaling hypothesis for domain walls, which has been put on a more rigorous 
footing in p5| , 27| . The scaling hypothesis, and the theory of dynamic scaling, can be also 



be applied to other extended topological defects such as cosmic strings. More generally, one 
can define a scaling exponent b, such that A = a/t b , and one of the goals of this paper is to 
measure both the exponent b and the amplitude a as accurately as possible. 

The first numerical simulations of this system were performed by Press, Ryden and 



Spergel |TT|, [17j who noted that in comoving coordinates the width of the wall shrinks as 
.R -1 , and so any numerical simulation with a lattice spacing fixed in comoving coordinates 
runs the risk of failing to resolve the domain wall at late times. They showed that ignoring 
the R dependence of the /i 2 parameter did not substantially affect the dynamics of the walls, 
and we adopt the same approach. We also include a damping term to simulate cooling in 
the early stages of the evolution. Written in first order form, the field equations become 

= 7T (6) 

V 2 - A0 {f - /i 2 ) -(^ + v) (7) 



7T 



The equations may be rescaled to \i = 1, A = 1, so that the width of the wall is 1. These 
classical equations can be thought of as effective equations representing the long- wavelength 
dynamics of the quantum field when the occupation number is high. 

The results of JTT] , |TT| seemed to show that a good fit for the area scaling law was given 

by 

A<x\n(t)/t, (8) 

in both 2 and 3 dimensions, which is not predicted by the theory of dynamic scaling for 
domain walls p5| , |57fl . It is therefore important to check these numerical simulations, and 
after a decade of development in computer technology one can do much larger simulations 
in order to eliminate transient effects. Indeed, our largest simulation is performed on 3D 
grid of 4096 3 , which gives us a dynamic range of roughly three orders of magnitude between 
t and M _1 , the characteristic response time of the field |3B . 



III. NUMERICAL SIMULATIONS 



A. Evolution Algorithm 

We used 2 and 3D cubic lattices with a 2nd order discretization of the Laplacian operator. 
The evolution of the discretized system was effected with the Verlet or leapfrog algorithm, 
common in molecular dynamics and offering a simple but effective way of discretizing the 
equations of motion (|7|). The scheme is 

VT n+ i = 7T n _l + / (t, <j> n , 7Tn) At 
0n+l = 0n + TTn+iAt, 

where f (t,(fi n ,ii n ) is the right hand side of Eq. (0), and 7r n = (7r n+ i —n n i)/2. As / is linear 
in 7r n , the equation can easily be rearranged so that vr n+ i is on the left hand side. The 
simulations on which the data is based all have Ax = 0.3 and At = 0.1. 

We chose periodic boundary conditions, and restricted the length of the simulation to a 
maximum time equal to T to t = NAx/2At, the time required for two signals emitted from the 
same point and travelling in opposite directions to interfere with each other. This represents 
the time it takes for the field to "notice" the finite dimension of the lattice it resides in. 
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FIG. 1: Typical evolution of \(j>\ d 



B. Initial Conditions 



The main objective of the simulations was not to see the formation of the network, 
but the evolution at later times, and so it is not necessary to start with a proper thermal 
distribution. Instead, we try to imitate the configuration of the field in the potential at some 
temperature above the critical T > Tc, with a Gaussian distribution around zero in fc-space. 
We created a Gaussian distribution by adding 12 uniformly distributed random numbers 
of unit variance, and algorithm that is quick and easy to parallelize. By setting <pk=o = 
we set the average of the distribution to zero and the proceed with transforming the field 
to the ^-representation. This is an alternative to simulating the actual phase transition by 
deforming the symmetry-breaking potential. The Gaussian distribution was chosen to give 
a pointwise spatial variance of about 10 _1 yU, and the value of the volume-averaged field was 
always less than 10~ 8 . Bigger ranges for the Gaussian distribution were found to lead to 
instabilities unless high dissipation was used. This is due to the occasional large values of 
the field, which have high energy densities due to the quartic potential. After initialisation 
the field is left to roll slowly to the two minima and start oscillating in them. A typical 
evolution of \<f)\ 2 can be seen in Fig. [|. 

If the average field were not close to zero in the initial conditions distribute we would 
have found that the biased initial conditions would give a different time dependence in the 
area density [12, 13fl . Indeed, if one phase is selected in the initial distribution then the 
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evolution of the domain wall area density can be shown to follow the relation 



A oc -e- ct , 



(10) 



where c a positive real number. Percolation theory predicts that for a phase occupying a 
fraction less than a critical value p c , the infinite wall disappears. The dominant phase takes 
over at a characteristic time scale, and quickly fills the entire simulation volume. In our case 
p ~ 0.5 ± 10~ 7 , far above the percolation threshold which for a cubic lattice can be shown 
to be p c = 0.31. 

C. Beginning the Evolution 

In the beginning of the simulation dissipation is imposed in order for the field to sink in 
a controlled manner into the minima and to remove spurious high frequency modes. This 
is controlled by the parameter 77 in (0). The dissipation aids the formation of the domain 
wall network and as soon as this is formed, dissipation is turned off and the system is left to 
evolve freely. A small amount of time, roughly the system's characteristic time, is required 
for the field to adjust itself into the new equation with no dissipation and continue to evolve. 
For that reason any regression on data starts a few timesteps after the end of the dissipated 
period. Turning the dissipation off slowly seems to help the system adjust more quickly to 
the new conditions and smooths out the effects of this 'adjustment' period. More specifically, 
the end of the dissipation period T D and the beginning of the regression T R are specified in 
the program and their difference Tr — To is calculated. As soon as the system reaches To, 
dissipation at timestep n rj n is decreased till zero using a Gaussian like function 



where n is the timestep, rji the initial dissipation and A a real positive number controlling 
in more accuracy the length of this period. 

Dissipations ranged from 0.1 to 0.3 in the simulations and affected the system for roughly 
10%-30% of the total simulation time T tot , with regression on the results starting at 20%-40% 
of the total time and the coefficient A was taken to be in the range 1 to 1.5 . 




(11) 
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D. Wall Area Density Calculation 



The main purpose of the simulations was to check the power law for the evolution of 
the wall area density. A crude way of calculating the total wall area is counting the places 
where adjacent lattice points have opposite signs and multiplying by Ax or Ax 2 which is a 
rough estimation of the wall area at the "link" . More precisely, one should find an accurate 
discretization of the continuum area density operator (|5]). 

The calculation is made by finding two adjacent lattice points where the field has opposite 
signs and calculating the gradient of the field at the "link". For two such adjacent points 
{i,j, k} and {i — l,j, k} the gradient is calculated as follows. The numerical approximation 
to the gradient at the direction of the link is just 

4>ij,k ~ 4>i-l,j,k 



A; 



Ax 



whereas the gradients at the remaining two directions are taken by averaging over the gra- 
dients of nearby links; the j -component of the gradient for example would be 

A (h — — ( ^-i.J+l.fc ~ 0i— l,J-l,fc 4>i,j+l,k ~ 4>i,j-l,k \ 

j<P ~ 2 V 2Ax" 2Ax" ) ' 

One needs to account for the orientation of the area element at each link, otherwise one will 
end up overestimating the area |L7|, |29| . For domain walls in 3D the problem can be resolved 



by simply multiplying the lattice area estimate by a factor | fl2"5 |, while a similar argument 
to that given in ]29|] gives the factor 7r/4 in 2D: 



A = A Lat ^ (3D) (12) 
A = A L J- (2D) (13) 

The Figures show A Lat , while the Tables show A. 

IV. RESULTS 

A. Wall area density 

In all simulations the wall area data started to be taken shortly after the dissipation had 
been stopped. Fig. |2] shows the results from a 2D run, fitted to a power law 
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FIG. 2: Evolution of the wall area density A^ a t for a 2D simulation, N = 1024, Ax = 0.3, At = 0.1. 
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at 



-b 



(14) 



by regression. Five simulations for the same parameters were run and an average was 
taken, with the fit taken on the averaged area. The results for 2D and 3D simulations 
are presented in Tables |, U for Minkowski, radiation and matter-dominated Friedmann- 



Background 


Length scaling law 


Minkowski 

Radiation 

Matter 


0.765(0.227).r - 987 (°- 032 ) 
0.928(0.165).r°- 996 ( a018 ) 
1.145(0.227).r°- 992 ( a014 ) 



TABLE I: Area scaling laws for Minkowski, radiation-dominated, and matter-dominated FRW 
backgrounds in 2 dimensions. The results are derived from averaging over 5 simulations for each 
case, with Ax = 0.3, At = 0.1 and a 1024 2 lattice. 



Robert son- Walker backgrounds respectively. The biggest 3D simulation (N = 4096, Ax = 
0.3, At = 0.08) gave a = 0.980(±0.017), 6 = 0.985(±0.003). 

There is a suggestion from the simulations that the wall area decreases slightly slower 
than 6=1. However, this deviation from the 6=1 scaling could not be attributed to a 



relation of the form given in Eq. [g, as it has been suggested [II]. Plotting exp(ALat) x t 
against time shows that there seems to be no logarithmic term in the wall area evolution, 

Fig. a 



Background 


Area scaling law 


Minkowski 

Radiation 

Matter 


0.883(0.141).r°- 995 ( a026 ) 
0.925(0.125).t-°- 994 ( a013 ) 
0.963(0.122).r°- 997 ( a012 ) 



TABLE II: Area scaling laws for Minkowski, radiation-dominated, and matter-dominated FRW 
backgrounds in 3 dimensions. The results are derived from averaging over 5 simulations for each 
case, with Ax = 0.3, At = 0.1 and a 512 3 lattice. 
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FIG. 3: exp(ALatt) against time for the simulation of Fig. ||. A decay law Ai at ~ log(i)/i would 
show logarithmic increase on this graph. 



Both the power law and the coefficient a of Eq. pj]. present challenges to the analytic 



method for the calculation of the wall area density of Ref . |25| , which are compared in detail 
in 



27j . The power law is very close to the predicted value of b = 1, with good precision, but 



the coefficient a showed larger fluctuations between runs. 



V. CONCLUSIONS 

The numerical integration of the equations of motion for the the 4 model has given an 
insight to the dynamics involved in domain wall networks and has provided an accurate way 
to support the scaling solution predicted by theoretical computations. A power law with 
exponent very close to 1 was found to be the best solution according to the simulation data 
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with no evidence for a logarithmic term suggested in previous studies |TT] . 

The fact that a domain wall network shows this dynamic scaling over approximately three 
orders of magnitude in the parameter M£, the ratio between the local curvature radius and 
the wall width, is remarkable, and has important implications for cosmic string networks. 
Firstly, it is clear from the visualisations in Fig. |] (see also P0|) that the energy in the 
domain walls is very quickly transferred into propagating modes of the field: the formation 
and collapse of closed loops (2D) or surfaces (3D), expected in the standard picture of the 
evolution of wall networks ||, is rare. Indeed, in 2D, self-intersections to form closed loops 
must be very rare; if two segments of wall approach each other they must be generically 
curved away from the point of closest approach, and therefore the acceleration is in the 
direction which would tend to increase the separation. Nonetheless, 2D walls scale perfectly 
well, so it seems plausible in that case that energy is being transferred directly into radiation. 
If the amplitude of the oscillations is large enough they could appear to form tiny loops or 



"protoloops" M, 16 



We believe that our results add weight to the contention, first put forward in [IJJ], that 
extended defects (including cosmic strings in 3D) have a non-perturbative channel into 
massive radiation. At first sight this is difficult to square with the standard picture, in 
which walls and strings obey the Nambu-Goto equations of motion for large curvature radii, 
for in that case the total energy locked up in the extended defects is conserved, in the absence 
of an general relativistic effects such as an expanding background or gravitational radiation. 
It is certainly true that it is possible to find string trajectories which are very close to being 
solutions of the Nambu-Goto equations [15, 3~lj : however, the initial conditions have to be 
carefully prepared, and the existence of these trajectories does not preclude the existence 
of a non-perturbative radiative process for defect networks. Indeed, we maintain that our 
results are good evidence that there must be such a process. 
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FIG. 4: Six snapshots from a numerical simulation using the algorithm described in Sec. Ill, with 
time increasing left to right and top to bottom. The solid isosurfaces are surfaces of constant <f>, 
with the red representing small positive values and the blue small negative values. Semi-transparent 
isosurfaces are surfaces of constant momentum density. 
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APPENDIX A: LATTICE CORRECTION FACTOR FOR 2D WALL LENGTH 
DENSITY 



In Ref. [f2^] it was shown that the naive lattice estimate of the area density of random 
domain walls in 3D overcounts the continuum value by a factor of 3/2 in the limit that the 
radii of curvature are large compared with the lattice spacing. In this section we perform 
the analagous calculation for two dimensions, finding it to be 4/tt. 

The naive estimate is obtained by summing the length of all links containing a wall and 
dividing by the total volume. A link crossing a wall is defined to be one for which the values 
of the field on the sites at either end have opposite signs. One can immediately see this 
will overestimate the length, as one is approximating a smooth curve by a sequence of line 
segments parallel to the lattice vectors i and j. 

In the continuum the centre of a domain wall is described by a 2-dimensional curve r(cr). 
Consider a segment of length /, with I <C £, where £ is the correlation length of the curve. 
Writing Ar = r(Z) — r(0) = xi + yj, we see that the lattice approximation to the length is 

k a t = \x\ + \y\. (Al) 

In the limit that £/£ — ► 0, the continuum value of the length is I = \/(x 2 + y 2 ). Hence 

/ Lat = ^ = (|cosa| + | sin/31)/ (A2) 

where a and (3 = n/2 — a are the direction cosines of the vector Ar. The ratio between the 
lattice estimate of the length and the true length is obtained by averaging over all possible 
orientations of the lattice relative to Ar: 



\ — — [ da(\ cosa| + | sinal) = — , (A3) 

I I ZTT JO 7T 



as advertised. 

This calculation can be extended to arbitrary I through a more involved argument. Let 
us first define the correlation function 

C(r) = (0(r)0(O)>, (A4) 
which is assumed to be smooth at r = 0, so that 

C(r) = C(0) - \c'\Qy. (A5) 
13 



In Refs. [p4, 25| it is shown that the length density A of the locus of zeroes of a Gaussian 



random field in 2 dimensions is given by 

A 1 



C "<°> (A6) 



2^ C(0)' 

This is the value of the length density in the continuum. 

On the lattice, we must consider the probability that the values of the field at opposite 
ends of a link have opposite signs, in which case we can say that the link is occupied by a 
segment of domain wall. Let us call this probability p occ , which is 

p occ = 2P(>(x) > AND 0(x + iAx) < 0), (A7) 

where the factor of 2 accounts for the opposite case 0(x) > AND 0(x + iAx) < 0. The 
lattice estimate of the length density is then 

*- = If- < A8 > 

where the lattice spacing is Ax, and the factor 2 comes from the fact that there are twice 
as many links as sites in 2 dimensions. 
Suppose we now define 

P 12 (C(r 12 ), V) = P(0( Xl ) > V AND 0(x 2 ) < V), (A9) 

where V is an arbitrary threshold and ru = |xi — x 2 |. One can then almost trivially write 

It can be shown |32] that 

= 1 eX p ( V% \ (All) 

where Ci 2 = C(r 12 ). Hence the lattice estimate of the area density is 

Providing Ax is much smaller than the correlation length £, defined by £ 2 = \C(0)/C"(0)\, 
we see that 

2 



^Lat 

7T 



(A13) 



\ C(Q) ' 
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and hence, from Eq. (|A6|) , that A^ at = AA/ir. 
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